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The ionization fraction of cosmic hydrogen, left over from the big bang, 
provides crucial fossil evidence for when the first stars and quasar black holes 
formed in the infant universe. Spectra of the two most distant quasars known^ 
show nearly complete absorption of photons with wavelengths shorter than the 
Lya transition of neutral hydrogen, indicating that hydrogen in the intergalactic 
medium (IGM) had not been completely ionized at a redshift z ~ 6.3, about 
a billion years after the big bang. Here we show that the radii of influence 
of ionizing radiation from these quasars imply that the surrounding IGM had 
a neutral hydrogen fraction of tens of percent prior to the quasar activity, 
much higher than previous lower limits ^'^ of ~ 0.1%. When combined with 
the recent inference of a large cumulative optical depth to electron scattering 
after cosmological recombination from the WMAP data'^, our result suggests 
the existence of a second peak in the mean ionization history, potentially due 
to an early formation episode of the first stars. 

The detection of a Lya emitting galaxy at 2; ~ 6.56 led to the claim that the IGM 
must already be highly ionized at that redshift^. However, it was later shown^ that the 
presence of a small ionized region around the galaxy and broadening of the Lya line. 
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would allow detection of the line feature in a fully neutral IGM. Luminous quasars offer an 
alternative probe of the IGM and have higher luminosities, allowing acquisition of higher 
quality spectra. The spectra^'^ of SDSS J1030+0524 at ^ = 6.28 and SDSS J1148+5251 at 
z — 6.41 show transmitted flux down to wavelengths corresponding to the Lya resonance 
of hydrogen at z ~ 6.20 and z ~ 6.32 respectively, even though no flux is detected at 
somewhat greater distances from the quasars''. These redshift displacements correspond to 
velocity shifts of order ~ 3300 km s~^, and imply ionized regions surrounding the quasars 
with observed physical radii of i?obs = 4.5 Mpc and 4.7 Mpc, respectively. The inference is 
robust since high ionization broad-lines typically have velocity offsets from the true quasar 
redshift^ of only ^ 1000 km s~^. 

We model the evolution of a fully ionized region (the so-called Stromgren sphere) 
around a quasar that emits UV photons into a partially neutral IGM. Prior to the overlap 
phase of ionized hydrogen (H II ) regions that signals the end of the reionization epoch^, 
galaxies generated their own small ionized regions flUing a fraction of the IGM volume 
around the quasar. We define the neutral fraction, xhi, to be the mean fraction of hydrogen 
atoms (in a representative volume of the low-density IGM) that remain neutral prior to the 
quasar activity, and assume xm outside the Stromgren sphere to remain constant during its 
expansion. Since the luminosity of a bright quasar is much larger than the stellar luminosity 
of neighboring galaxies, a fully ionized quasar bubble expands into a partially ionized IGM 
composed of isolated H II regions. 

Even for the deepest available spectra, with lower limits on the Lya ( Gunn-Peterson}'^) 
optical depth of tqp ^ 26, the IGM need only have a neutral fraction of ~ 10~^ 
in order to produce the observed absorption blueward of Lya^. However, the value 
of i?obs ~ 4.5Mpc may be used to constrain larger values of the neutral fraction in 
the range 0.001 < xhi < 1.0. In the early expansion phase, the dependence of the 
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physical radius of the Stromgren sphere, Rp, on the quasar age^'^^, ^age; ma-y be crudely 
approximated as Rp ~ 7x^1^^ [i^agc/lO^yrj^'^^Mpc, given the production rate of ionizing 
photons that is characteristic of SDSS J1148+5251 and SDSS J1030+0524. Clearly, 
the Stromgren sphere grows larger when embedded in an IGM with a small neutral 
fraction. Defining i?max to be the maximum radius achieved in a fully neutral IGM 
by the end of the quasar lifetime, tit, the resulting a priori probability distribution is, 
dP/dR,\^^, = 3xHii?^/i?La. for Rp < R^^x^l^K In addition, the neutral fraction must 
be smaller than a^max = min (^[tit/lO^ears] [i?obs/7Mpc]~^ , to allow the growth of the 
sphere to i?obs within tu. Using a flat logarithmic prior for xm, we flnd the cumulative 
a posteriori probability distribution, P{< xm) = min(a;Hi, a;niax)/^max- Since estimates of 
quasar hfetimes^^ are bracketed in the range 10^ — 10^ years, we get Xmax ~ 1 and conclude 
that Xm ;^ 0.1 (0.01) with 90% (99%) confidence. 

To better quantify this simple argument, we write the general equation for the 
relativistic expansion of the co-moving radius [R— (1 + z)Rp] of the quasar H II region^ in 
an IGM with a neutral filling fraction xm (fixed by other ionizing sources), 

dR_ i F^Nion - o^BCF^xni (n^)'' (1 + zf (f i^^) 

^ F^iVion + 47ri?2 (1 + z) cF^xmn^ 

where c is the speed of light, is the mean number density of protons at ^ = 0, 
ckb = 2.6 X 10~^^cm^s~^ is the case-B recombination coefficient at the characteristic 
temperature of lO^K, and A^ion is the rate of ionizing photons crossing a shell at the radius 
of the H II region at time t. We use the distribution derived from numerical simulations 
for the over-densities in gas clumps^^, and calculate the mean free path o?(Ac) for ionizing 
photons^^'^^. Following Barkana & Loeb^^, we then find the critical overdensity (Ac) at 
which a fraction = 0.5 = exp [— i?p/(i(Ac)], of photons emitted do not encounter an 
overdensity larger than Ac within the H II region. We also compute the mass fraction Fm 
(~ 1) of gas within Rp that is at over-densities lower than Ac. Finally, we calculate the 
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clumping factor in the ionized regions, C{R) = (A^)/(A)^, where the angular brackets 
denote an average over all regions with A < Ap. For xhi = 1 and dRp/dt <^ c, equation (1) 
reduces to its well-known form^^'^^~^^. The expansion of the quasar Stromgren sphere 
would change in an overdense region of the IGM. However the density contrast due to infall 
has been shown-*^^ to be small (~ 1%) on scales of several Mpc around the massive hosts 
of the bright SDSS quasars at 2; ~ 6. Numerical experiments showed our results to be 
insensitive to infall. The emission rate of ionizing photons A^ion in equation (1) is computed 
at f — t — ideiay to accouut for the finite light travel time between the source and the 
ionization front. The delay tdeiay is derived from the relation R = It^-t^^^^^ cdt' [1 + z{t')], 
where is the time when the photon crosses the co-moving radius R. The use of the 
Telfer et al.^^ spectral template imphes Ni^^ = 6.5 x lO^^s-^ for SDSS J1030+0524 and 
TVion = 10.0 X lO^^s-^ for SDSS J1148+5251, while the template of Elvis et al.^'^ imphes 
lower values of iVion ~ 1-3 x 10^''s~^ and iVjon ~ 2.0 x 10^'^s"^, respectively^ 

Our model for the evolution of the Stromgren sphere includes the quasar formation 
history. We compute the time dependent ionizing flux governed by black-hole growth 
through accretion and mergers. Following the observational inference that the relation 
between bulge velocity dispersion and black hole mass Mbh does not evolve with redshift^°, 
we extrapolate the present-day Mbh-Mhaio relation^-^ using the dependence of virial 
velocity on redshift^ to obtain Mhaio = 1-5 x IO^^^^^q (Mbh/lO^Me)^/^ [(1 + z)/7]~^^^. The 
luminosities of the SDSS quasars imply black- hole masses of ~ 2 x lO^M© accreting at their 
Eddington rate^, leading to inferred host dark-matter halos of Mhaio ~ 2 x lO^^M©. The 
inference of such a massive halo at this early epoch is supported by the signature of gas 
infall in the spectra of high redshift quasars^'', and the inference of a large molecular mass 
and velocity width in the host galaxy of SDSS J1148-|-5251^. 

We begin with two dark-matter halos of comparable mass Mi and M2 that merge to 
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form the host dark-matter halo with a mass Mhaio = (-^i + -^2) at a redshift corresponding 
to one quasar hfetime prior the observed redshift, and construct merger trees for each 
of the progenitor halos using the algorithm described by Volonteri et al.^^. We use the 
Mbh-Mhaio relation to find the mass of the black-hole at the center of each dark-matter halo 
in the merger tree. When there is a major merger (having a progenitor mass ratio < 2), we 
assume that the black-holes merge and that the accretion of the mass deficit relative to the 
Mbh-Mhaio relation produces the limiting Eddington luminosity with a radiative efficiency of 
e = 10%. The resulting quasar hfetime is tit = 4 x 10^(e/0.1) In [mI/^^J{mI'^ + M^'^)] years. 
The ionizing photon emission rate N{t) is taken to have an exponential time dependence, 
-^ion(t) = A/'oe"*/*'*, with Nq oc Mbh- We assume that all quasar episodes in the merger 
tree contribute to the volume of the observed Stromgren sphere, which is centered on the 
most massive halo in the tree. A similar prescription for quasar activity was shown to be 
consistent with the observed number counts of high redshift quasars^^. For a major merger 
resulting in the observed quasar activity, the lifetime given by this approach is ~ 10^ years. 
This compares favorably with a variety observational estimates of quasar lifetimes^^"^^. To 
cover the full range of uncertainty we consider modifications to our fiducial model in which 
the quasar lifetime is multiplied by a factor between 0.1-10, resulting in lifetimes for the 
observed quasars of /it^it = 10^-10^ years. 

We have produced 300 realizations of the merger tree, and computed the evolution of 
the Stromgren sphere in each case for the full range of neutral fractions 0.001 < xhi ^ 1.0, 
allowed by the Gunn- Peterson optical depth^ . We find that quasar activity associated with 
the hierarchical build up of the host galaxy produces an H II region with a typical size of 

— 1/3 

~ 2xjjj Mpc, comparable to the observed radii if xhi is of order unity. Figure 1 shows 
the conditional a priori probability distributions of for quasars like SDSS J1030-I-0524, 
assuming xhi = 1 and two different values of /it. The observed radii are consistent with 
lifetimes close to the fiducial case of /it = 1. The plotted distribution of implies that the 
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hierarchical evolution of early quasars leads to a dearth of small Stromgren sphere radii for 
the observed SDSS quasars. This is in contrast to the simplified monolithic model of quasar 
formation^^'^, for which the H II region has zero volume at the time when the quasar turns 
on [R{ta,ge = 0) = 0] and the distribution extends down to i?p = 0. 

To constrain the neutral fraction, we have computed the likelihood of observing 
Rohs = 4.5Mpc around SDSS J1030+0524 and R^hs = 4.7Mpc around SDSS J1148+5251 
as a function of xhi and f^. The results are plotted in Figure 2. The upper panels show 
the locus of most likely values (thick line) as well as likelihood contours at 0.1 of the peak 
value (dashed lines). The fiducial lifetime (/it) favors xm ~ 1, while with /it -C 1 the 
distributions for Rp lie substantially below the observed value of 4.5Mpc, making smaller 
values of xm more likely. Extrapolation of the most likely contour for the Elvis et al.^^ 
template yields /it ~ 2a;Hi. This implies that xm ~ 10~^ would require a lifetime as short 
as 2 X 10^ years, which is ruled out by variabihty properties of quasars in SDSS^^. The 
a posteriori probability distributions for xm given /it are plotted in the lower panels of 
Figure 2 and robustly yield the constraint xhi ^ 0.01. For /it ^ 0.3, we find Xm ^0.1 and 
Xm k, 0.4 at the 90% level for the Elvis et al.^^ and Telfer et al.^® spectra, respectively. The 
fiducial model with /it = 1 yields corresponding constraints of xhi ^ 0.3 and xhi ^ 0.6. 

The inference of a large neutral fraction at z ~ 6.3 presents a challenge to theories 
of cosmological reionization when combined with the large optical depth^ to electron 
scattering after cosmological recombination, Tes = 0.17 ± 0.04. Consider a toy model in 
which the universe was partially reionized at a high redshift 2;reion,i, leaving a neutral 
fraction xm until complete reionization was reached at -Zreion,!! ~ 6.25. The optical depth 
is then = 0.04 + 0.002(1 - xni) [(1 + -^reionj)^/^ - 19.5] . If the IGM ionization fraction 
increased monotonically, then given Tes > 0.13, the universe needed to be reionized earlier 
than 2;reion,i ~ 18 or Zreion,! ~ 24 assumiug the Elvis et al.^^ and Telfer et al.^* spectra. 
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respectively. If r^s — 0.17 and x-fji > 0.7, as implied by the Telfer et al.^^ spectrum, then 
a monotonic reionization history requires significant reionization at an implausibly high 
redshift {z ^ 30). In this case, a more plausible alternative is a non- monotonic history 
with an early reionization peak, possibly due to the formation of massive population-Ill 
stars^''''^^. 

Correspondence and requests for materials to Abraham Loeb. 
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Fig. 1. — Predicted probability for observing different radii of the ionized region around tlie 
quasar SDSS J1030+0524. The differential a priori probability distribution is plotted for 
xm — 1 and two quasar lifetimes: /it = 0.1 (dashed) and /it = 1.0 (solid). The light and 
dark curves show results calculated based on the Elvis et al.^^ and Teller et al.^^ template 
spectra, respectively. The distributions may also be used to approximate the behavior with 
xhi < 1 by replacing Rp with RpX^^^. Values of Rp significantly larger than observed still 
allow transmission of Lya flux at a detectable level because: (i) the optical depth due to the 
damping wing of the IGM^^ is only important near the boundary of the H II region, and (ii) 
the optical depth at Rp due to resonant absorption within the H II region^^ is in the range 
Treson = (l ~8) X (-Rp/4.5 Mpc)^. A deuse H I cloud near Robs could produce a large damping 
wing and lead to a systematic underestimation of Rp. However, such a cloud would need a 
column density ^ 10^^ cm~^ (;:^ 10^^ cm~^) to produce a Lya optical depth > 26 over 10% 
(30%) of the spectral range covered by the H II region. This requires a damped Lya absorber, 
whose existence is highly improbable within the narrow redshift interval under consideration, 
Az ~ 0.1. Our calculation assumes that the ionization front is thin. The spectrum-averaged 
mean-free path for the ionizing quasar photons is A ~ 1.5x^1(1 + 2;)~^Mpc, which may 
be compared to the bubble radius Rp ~ A.Bx^y^Mpc, to yield the fractional thickness of 
the ionization front (\/Rp) ~ lO'^ (i?p/4.5Mpc)"^ ^hi^^ [(1 + z)/7.3]"'l We also note that 
the ionized region need not be spherical; if the quasar radiation is beamed, its luminosity 
per unit solid- angle along the observer's line-of-sight is still measured, and equation (1) 
therefore describes the extent of the ionized region along the line-of-sight. Throughout our 
calculations, we have adopted the best-fit cosmological parameters derived from the WMAP 
data^°. 
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Fig. 2. — Likelihood for the inferred neutral fraction of the IGM, assuming different 
quasar lifetimes. We show contours of likelihood, L, for xm and /it {top) and a posteriori 
cumulative probability distributions for xm (bottom). Cumulative distributions are shown 
for two different quasar lifetimes: /it = 0.1 (dashed) and /it = 1.0 (solid). The light 
and dark lines show results where the quasar ionizing photon rate was specified by the 
Elvis et al.^^ and Telfer et al.^® spectra, respectively. The a posteriori probability is 
oc 

the prior probability for xhi, assumed to be fiat in logarithmic bins within the range 
0.001 ^ xhi ^ 1-0, and and are the a priori probability distributions for Rp 

in quasars like SDSS J1030+0524 and SDSS J1148+5251, respectively A flat logarithmic 
prior is the natural choice for because x^i may be thought of as a ratio of independent 

quantities (the number of ionizing photons and the number of baryons) that would themselves 
have linearly distributed prior probabilities. The alternative use of a linear prior in xm. leads 
to more stringent limits than those presented here. 
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